# Nighttime lights and land cover data are in a single chunk to avoid external pointer errors when rendering
#---
# Nighttime lights
#---
# Subset nighttime lights rasters to each city
city_names <- c("flagstaff", "homerglen", "beecave", "scranton",
"belvidere", "clanton")
cities <- list(flagstaff, homerglen, beecave, scranton, belvidere, clanton)
city_rasters <- vector("list", length = length(cities))
# Load files for each year
for (i in 1:length(years)){
pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",
years[i], ".tif", sep = ""), repo = "cstaebell/final-project-cstaebell")
rast_file <- paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",
years[i], ".tif", sep = "")
year_rast <- rast(rast_file)
# Loop over cities
for (j in 1:length(cities)) {
# Crop and mask to city extent
current_rast <- year_rast |>
terra::crop(cities[[j]]) |>
mask(cities[[j]])
# Create list entry during first iteration, then append in subsequent iterations
if (i == 1) {
city_rasters[[j]] <- list(current_rast)
} else {
city_rasters[[j]][[length(city_rasters[[j]]) + 1]] <- current_rast
}
}
}
# Stack city lists
city_rasters <- lapply(city_rasters, rast)
# Create separate rasters for each city
city_var_names <- paste0(city_names, "_rast")
for (c in 1:length(city_var_names)) {
assign(city_var_names[c], city_rasters[[c]])
}
# Calculate annual means
annual_means <- vector("list", length = length(cities))
for (i in 1:length(city_rasters)) {
ntl_vals <- data.frame(values(city_rasters[[i]])) |>
summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE)))
annual_means[[i]] <- ntl_vals
}
# Create mean dataframes for each city
city_df_var_names <- paste0(city_names, "_means")
for (i in 1:length(annual_means)) {
df <- annual_means[[i]] |>
pivot_longer(cols = t2014:t2024) |>
mutate(name = years) |>
rename(year = name, mean_ntl = value)
assign(city_df_var_names[i], df)
}
#---
# Land cover data
#---
# Download land cover data
pb_download("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif",
repo = "cstaebell/final-project-cstaebell")
pb_download("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif",
repo = "cstaebell/final-project-cstaebell")
land_use_rast_2019 <- rast("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif")
land_use_rast_2020 <- rast("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif")
# Subset land cover rasters to each city
land_use_rasters <- vector("list", length = length(cities))
urban_ntl_dfs <- vector("list", length = length(cities))
urban_ntl_vals <- vector("list", length = length(cities))
group_assign <- c("Dark Sky", "Dark Sky", "Dark Sky", "Comparison", "Comparison", "Comparison")
pair_assign <- c("A", "B", "C", "A", "B", "C")
for (i in 1:length(cities)) {
# Specify which data to use (2019 for all cities but Clanton)
if (i == 6) {
land_use_rast <- land_use_rast_2020
} else {
land_use_rast <- land_use_rast_2019
}
# Crop and mask to city extent
current_rast <- land_use_rast |>
terra::crop(cities[[i]]) |>
mask(cities[[i]])
# Save to list for future use
land_use_rasters[[i]] <- current_rast
# Mask all but urban land use (land use code = 13)
urban_area <- current_rast
urban_area[urban_area != 13] <- NA
city_urban <- city_rasters[[i]] |>
mask(urban_area)
# Calculate mean NTL for urban areas
urban_ntl_dfs[[i]] <- data.frame(values(city_urban)) |>
summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |>
pivot_longer(cols = t2014:t2024) |>
mutate(name = years) |>
rename(year = name, mean_urban_ntl = value)
# Compile values for statistical tests
# put these outside the loop when done. also replace paste0 with city_name
urban_ntl_vals[[i]] <- data.frame(values(city_urban)) |>
pivot_longer(cols = t2014:t2024) |>
filter(!is.na(value)) |>
mutate(city = city_names[[i]], group = group_assign[i], pair = pair_assign[i])
}
# Create separate dataframes from list object
urban_ntl_var_names <- paste0(city_names, "_urban_ntl")
for (i in 1:length(urban_ntl_dfs)) {
assign(urban_ntl_var_names[i], urban_ntl_dfs[[i]])
}
# Prepare data for sample land cover plots
# Define land cover encoding
Land_Cover_Type_1 <- c(
"Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2,
"Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Mixed forest" = 5,
"Closed shrublands" = 6, "Open shrublands" = 7, "Woody savannas" = 8, "Savannas" = 9,
"Grasslands" = 10, "Permanent wetlands" = 11, "Croplands" = 12, "Urban & built-up" = 13,
"Cropland/Natural vegetation mosaic" = 14, "Snow & ice" = 15, "Barren/Sparsely vegetated" = 16
)
lcd <- data.frame(
ID = Land_Cover_Type_1,
landcover = names(Land_Cover_Type_1),
stringsAsFactors = FALSE
)
# Create data frames for plotting
scranton_lulc <- as.data.frame(land_use_rasters[[4]], xy = TRUE) |>
left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID))
# Since Scranton has more landcover types, use its types to join with Flagstaff for a consistent color scheme
lcd_reduced <- lcd |>
filter(landcover %in% unique(scranton_lulc$landcover))
flagstaff_lulc <- as.data.frame(land_use_rasters[[1]], xy = TRUE) |>
right_join(lcd_reduced, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID))